// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_RANDOMSETTER_H
#define EIGEN_RANDOMSETTER_H

#if defined(EIGEN_GOOGLEHASH_SUPPORT)
// Ensure the ::google namespace exists, required for checking existence of
// ::google::dense_hash_map and ::google::sparse_hash_map.
namespace google {
}
#endif

namespace Eigen {

/** Represents a std::map
  *
  * \see RandomSetter
  */
template <typename Scalar> struct StdMapTraits
{
    typedef int KeyType;
    typedef std::map<KeyType, Scalar> Type;
    enum
    {
        IsSorted = 1
    };

    static void setInvalidKey(Type&, const KeyType&) {}
};

#ifdef EIGEN_UNORDERED_MAP_SUPPORT
/** Represents a std::unordered_map
  *
  * To use it you need to both define EIGEN_UNORDERED_MAP_SUPPORT and include the unordered_map header file
  * yourself making sure that unordered_map is defined in the std namespace.
  *
  * For instance, with current version of gcc you can either enable C++0x standard (-std=c++0x) or do:
  * \code
  * #include <tr1/unordered_map>
  * #define EIGEN_UNORDERED_MAP_SUPPORT
  * namespace std {
  *   using std::tr1::unordered_map;
  * }
  * \endcode
  *
  * \see RandomSetter
  */
template <typename Scalar> struct StdUnorderedMapTraits
{
    typedef int KeyType;
    typedef std::unordered_map<KeyType, Scalar> Type;
    enum
    {
        IsSorted = 0
    };

    static void setInvalidKey(Type&, const KeyType&) {}
};
#endif  // EIGEN_UNORDERED_MAP_SUPPORT

#if defined(EIGEN_GOOGLEHASH_SUPPORT)

namespace google {

    // Namespace work-around, since sometimes dense_hash_map and sparse_hash_map
    // are in the global namespace, and other times they are under ::google.
    using namespace ::google;

    template <typename KeyType, typename Scalar> struct DenseHashMap
    {
        typedef dense_hash_map<KeyType, Scalar> type;
    };

    template <typename KeyType, typename Scalar> struct SparseHashMap
    {
        typedef sparse_hash_map<KeyType, Scalar> type;
    };

}  // namespace google

/** Represents a google::dense_hash_map
  *
  * \see RandomSetter
  */
template <typename Scalar> struct GoogleDenseHashMapTraits
{
    typedef int KeyType;
    typedef typename google::DenseHashMap<KeyType, Scalar>::type Type;
    enum
    {
        IsSorted = 0
    };

    static void setInvalidKey(Type& map, const KeyType& k) { map.set_empty_key(k); }
};

/** Represents a google::sparse_hash_map
  *
  * \see RandomSetter
  */
template <typename Scalar> struct GoogleSparseHashMapTraits
{
    typedef int KeyType;
    typedef typename google::SparseHashMap<KeyType, Scalar>::type Type;
    enum
    {
        IsSorted = 0
    };

    static void setInvalidKey(Type&, const KeyType&) {}
};
#endif

/** \class RandomSetter
  *
  * \brief The RandomSetter is a wrapper object allowing to set/update a sparse matrix with random access
  *
  * \tparam SparseMatrixType the type of the sparse matrix we are updating
  * \tparam MapTraits a traits class representing the map implementation used for the temporary sparse storage.
  *                  Its default value depends on the system.
  * \tparam OuterPacketBits defines the number of rows (or columns) manage by a single map object
  *                        as a power of two exponent.
  *
  * This class temporarily represents a sparse matrix object using a generic map implementation allowing for
  * efficient random access. The conversion from the compressed representation to a hash_map object is performed
  * in the RandomSetter constructor, while the sparse matrix is updated back at destruction time. This strategy
  * suggest the use of nested blocks as in this example:
  *
  * \code
  * SparseMatrix<double> m(rows,cols);
  * {
  *   RandomSetter<SparseMatrix<double> > w(m);
  *   // don't use m but w instead with read/write random access to the coefficients:
  *   for(;;)
  *     w(rand(),rand()) = rand;
  * }
  * // when w is deleted, the data are copied back to m
  * // and m is ready to use.
  * \endcode
  *
  * Since hash_map objects are not fully sorted, representing a full matrix as a single hash_map would
  * involve a big and costly sort to update the compressed matrix back. To overcome this issue, a RandomSetter
  * use multiple hash_map, each representing 2^OuterPacketBits columns or rows according to the storage order.
  * To reach optimal performance, this value should be adjusted according to the average number of nonzeros
  * per rows/columns.
  *
  * The possible values for the template parameter MapTraits are:
  *  - \b StdMapTraits: corresponds to std::map. (does not perform very well)
  *  - \b GnuHashMapTraits: corresponds to __gnu_cxx::hash_map (available only with GCC)
  *  - \b GoogleDenseHashMapTraits: corresponds to google::dense_hash_map (best efficiency, reasonable memory consumption)
  *  - \b GoogleSparseHashMapTraits: corresponds to google::sparse_hash_map (best memory consumption, relatively good performance)
  *
  * The default map implementation depends on the availability, and the preferred order is:
  * GoogleSparseHashMapTraits, GnuHashMapTraits, and finally StdMapTraits.
  *
  * For performance and memory consumption reasons it is highly recommended to use one of
  * Google's hash_map implementations. To enable the support for them, you must define
  * EIGEN_GOOGLEHASH_SUPPORT. This will include both <google/dense_hash_map> and
  * <google/sparse_hash_map> for you.
  *
  * \see https://github.com/sparsehash/sparsehash
  */
template <typename SparseMatrixType,
          template <typename T> class MapTraits =
#if defined(EIGEN_GOOGLEHASH_SUPPORT)
              GoogleDenseHashMapTraits
#elif defined(_HASH_MAP)
              GnuHashMapTraits
#else
              StdMapTraits
#endif
          ,
          int OuterPacketBits = 6>
class RandomSetter
{
    typedef typename SparseMatrixType::Scalar Scalar;
    typedef typename SparseMatrixType::StorageIndex StorageIndex;

    struct ScalarWrapper
    {
        ScalarWrapper() : value(0) {}
        Scalar value;
    };
    typedef typename MapTraits<ScalarWrapper>::KeyType KeyType;
    typedef typename MapTraits<ScalarWrapper>::Type HashMapType;
    static const int OuterPacketMask = (1 << OuterPacketBits) - 1;
    enum
    {
        SwapStorage = 1 - MapTraits<ScalarWrapper>::IsSorted,
        TargetRowMajor = (SparseMatrixType::Flags & RowMajorBit) ? 1 : 0,
        SetterRowMajor = SwapStorage ? 1 - TargetRowMajor : TargetRowMajor
    };

public:
    /** Constructs a random setter object from the sparse matrix \a target
      *
      * Note that the initial value of \a target are imported. If you want to re-set
      * a sparse matrix from scratch, then you must set it to zero first using the
      * setZero() function.
      */
    inline RandomSetter(SparseMatrixType& target) : mp_target(&target)
    {
        const Index outerSize = SwapStorage ? target.innerSize() : target.outerSize();
        const Index innerSize = SwapStorage ? target.outerSize() : target.innerSize();
        m_outerPackets = outerSize >> OuterPacketBits;
        if (outerSize & OuterPacketMask)
            m_outerPackets += 1;
        m_hashmaps = new HashMapType[m_outerPackets];
        // compute number of bits needed to store inner indices
        Index aux = innerSize - 1;
        m_keyBitsOffset = 0;
        while (aux)
        {
            ++m_keyBitsOffset;
            aux = aux >> 1;
        }
        KeyType ik = (1 << (OuterPacketBits + m_keyBitsOffset));
        for (Index k = 0; k < m_outerPackets; ++k) MapTraits<ScalarWrapper>::setInvalidKey(m_hashmaps[k], ik);

        // insert current coeffs
        for (Index j = 0; j < mp_target->outerSize(); ++j)
            for (typename SparseMatrixType::InnerIterator it(*mp_target, j); it; ++it)
                (*this)(TargetRowMajor ? j : it.index(), TargetRowMajor ? it.index() : j) = it.value();
    }

    /** Destructor updating back the sparse matrix target */
    ~RandomSetter()
    {
        KeyType keyBitsMask = (1 << m_keyBitsOffset) - 1;
        if (!SwapStorage)  // also means the map is sorted
        {
            mp_target->setZero();
            mp_target->makeCompressed();
            mp_target->reserve(nonZeros());
            Index prevOuter = -1;
            for (Index k = 0; k < m_outerPackets; ++k)
            {
                const Index outerOffset = (1 << OuterPacketBits) * k;
                typename HashMapType::iterator end = m_hashmaps[k].end();
                for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it != end; ++it)
                {
                    const Index outer = (it->first >> m_keyBitsOffset) + outerOffset;
                    const Index inner = it->first & keyBitsMask;
                    if (prevOuter != outer)
                    {
                        for (Index j = prevOuter + 1; j <= outer; ++j) mp_target->startVec(j);
                        prevOuter = outer;
                    }
                    mp_target->insertBackByOuterInner(outer, inner) = it->second.value;
                }
            }
            mp_target->finalize();
        }
        else
        {
            VectorXi positions(mp_target->outerSize());
            positions.setZero();
            // pass 1
            for (Index k = 0; k < m_outerPackets; ++k)
            {
                typename HashMapType::iterator end = m_hashmaps[k].end();
                for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it != end; ++it)
                {
                    const Index outer = it->first & keyBitsMask;
                    ++positions[outer];
                }
            }
            // prefix sum
            StorageIndex count = 0;
            for (Index j = 0; j < mp_target->outerSize(); ++j)
            {
                StorageIndex tmp = positions[j];
                mp_target->outerIndexPtr()[j] = count;
                positions[j] = count;
                count += tmp;
            }
            mp_target->makeCompressed();
            mp_target->outerIndexPtr()[mp_target->outerSize()] = count;
            mp_target->resizeNonZeros(count);
            // pass 2
            for (Index k = 0; k < m_outerPackets; ++k)
            {
                const Index outerOffset = (1 << OuterPacketBits) * k;
                typename HashMapType::iterator end = m_hashmaps[k].end();
                for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it != end; ++it)
                {
                    const Index inner = (it->first >> m_keyBitsOffset) + outerOffset;
                    const Index outer = it->first & keyBitsMask;
                    // sorted insertion
                    // Note that we have to deal with at most 2^OuterPacketBits unsorted coefficients,
                    // moreover those 2^OuterPacketBits coeffs are likely to be sparse, an so only a
                    // small fraction of them have to be sorted, whence the following simple procedure:
                    Index posStart = mp_target->outerIndexPtr()[outer];
                    Index i = (positions[outer]++) - 1;
                    while ((i >= posStart) && (mp_target->innerIndexPtr()[i] > inner))
                    {
                        mp_target->valuePtr()[i + 1] = mp_target->valuePtr()[i];
                        mp_target->innerIndexPtr()[i + 1] = mp_target->innerIndexPtr()[i];
                        --i;
                    }
                    mp_target->innerIndexPtr()[i + 1] = internal::convert_index<StorageIndex>(inner);
                    mp_target->valuePtr()[i + 1] = it->second.value;
                }
            }
        }
        delete[] m_hashmaps;
    }

    /** \returns a reference to the coefficient at given coordinates \a row, \a col */
    Scalar& operator()(Index row, Index col)
    {
        const Index outer = SetterRowMajor ? row : col;
        const Index inner = SetterRowMajor ? col : row;
        const Index outerMajor = outer >> OuterPacketBits;  // index of the packet/map
        const Index outerMinor = outer & OuterPacketMask;   // index of the inner vector in the packet
        const KeyType key = internal::convert_index<KeyType>((outerMinor << m_keyBitsOffset) | inner);
        return m_hashmaps[outerMajor][key].value;
    }

    /** \returns the number of non zero coefficients
      *
      * \note According to the underlying map/hash_map implementation,
      * this function might be quite expensive.
      */
    Index nonZeros() const
    {
        Index nz = 0;
        for (Index k = 0; k < m_outerPackets; ++k) nz += static_cast<Index>(m_hashmaps[k].size());
        return nz;
    }

protected:
    HashMapType* m_hashmaps;
    SparseMatrixType* mp_target;
    Index m_outerPackets;
    unsigned char m_keyBitsOffset;
};

}  // end namespace Eigen

#endif  // EIGEN_RANDOMSETTER_H
